# =============================================================================
# COMPREHENSIVE BOUNDARY CONDITION TEST SUITE
# =============================================================================

using Test
using PDEJulia
using PDEJulia.Grids
using LinearAlgebra
using SparseArrays

@testset "Boundary Conditions" begin
    
    # =============================================================================
    # BASIC BOUNDARY CONDITION CREATION AND VALIDATION
    # =============================================================================
    
    @testset "Basic BC Creation" begin
        # Test all basic BC types
        @test PeriodicBC() isa PeriodicBC
        @test bc_type(PeriodicBC()) == :periodic
        
        @test DirichletBC(1.0) isa DirichletBC
        @test DirichletBC(1.0).value == 1.0
        @test DirichletBC(0.0).homogeneous == true
        @test DirichletBC(1.0).homogeneous == false
        @test bc_type(DirichletBC(1.0)) == :dirichlet
        
        @test NeumannBC(0.5) isa NeumannBC
        @test NeumannBC(0.5).value == 0.5
        @test NeumannBC().value == 0.0  # Default
        @test bc_type(NeumannBC()) == :neumann
        
        @test RobinBC(1.0, 0.1, 2.0) isa RobinBC
        @test RobinBC(1.0, 0.1).value == 0.0  # Default value
        @test bc_type(RobinBC(1.0, 0.1)) == :robin
        
        # Test Robin BC validation
        @test_throws ArgumentError RobinBC(0.0, 0.0)  # Both coefficients zero
        
        mixed_bc = MixedBC(DirichletBC(1.0), NeumannBC(0.0))
        @test mixed_bc isa MixedBC
        @test bc_type(mixed_bc) == :mixed
        
        # Test MixedBC validation
        @test_throws ArgumentError MixedBC(PeriodicBC(), DirichletBC(1.0))
    end
    
    @testset "Time-Dependent BC" begin
        # Test time-dependent BC wrapper
        td_bc = TimeDependentBC(t -> DirichletBC(sin(t)))
        @test td_bc isa TimeDependentBC
        @test bc_type(td_bc) == :time_dependent
        
        # Check that initial BC is set
        @test td_bc.current_bc[] isa DirichletBC
        @test td_bc.current_bc[].value ≈ sin(0.0)
    end
    
    # =============================================================================
    # BOUNDARY CONDITION PARSING
    # =============================================================================
    
    @testset "BC Parsing" begin
        # Test string parsing
        @test parse_boundary_condition("periodic") isa PeriodicBC
        
        dirichlet_bc = parse_boundary_condition("dirichlet:1.5")
        @test dirichlet_bc isa DirichletBC
        @test dirichlet_bc.value ≈ 1.5
        
        neumann_bc = parse_boundary_condition("neumann:-0.5")
        @test neumann_bc isa NeumannBC
        @test neumann_bc.value ≈ -0.5
        
        robin_bc = parse_boundary_condition("robin:1.0,0.1,2.0")
        @test robin_bc isa RobinBC
        @test robin_bc.alpha ≈ 1.0
        @test robin_bc.beta ≈ 0.1
        @test robin_bc.value ≈ 2.0
        
        # Test default value for Robin
        robin_bc_default = parse_boundary_condition("robin:1.0,0.1")
        @test robin_bc_default.value ≈ 0.0
        
        # Test error cases
        @test_throws ArgumentError parse_boundary_condition("unknown:1.0")
        @test_throws ArgumentError parse_boundary_condition("robin:1.0")  # Missing beta
        
        # Test vector parsing
        bc_strings = ["periodic", "dirichlet:1.0", "neumann:0.0"]
        bcs = parse_boundary_conditions(bc_strings)
        @test length(bcs) == 3
        @test bcs[1] isa PeriodicBC
        @test bcs[2] isa DirichletBC
        @test bcs[3] isa NeumannBC
    end
    
    # =============================================================================
    # GRID INTEGRATION TESTS
    # =============================================================================
    
    @testset "Grid BC Integration" begin
        # Create test grids
        grid_1d = UnitGrid(32, (0.0, 1.0))
        grid_2d = CartesianGrid([32, 16], [(0.0, 1.0), (0.0, 0.5)])
        
        # Test setting BCs
        bcs_1d = [NeumannBC(0.0)]
        set_boundary_conditions!(grid_1d, bcs_1d)
        @test get_boundary_conditions(grid_1d) == bcs_1d
        
        bcs_2d = [PeriodicBC(), DirichletBC(1.0)]
        set_boundary_conditions!(grid_2d, bcs_2d)
        @test get_boundary_conditions(grid_2d) == bcs_2d
        
        # Test dimension mismatch
        @test_throws ArgumentError set_boundary_conditions!(grid_1d, [NeumannBC(), DirichletBC()])\n        \n        # Test individual axis access\n        @test get_boundary_condition(grid_2d, 1) isa PeriodicBC\n        @test get_boundary_condition(grid_2d, 2) isa DirichletBC\n        @test_throws BoundsError get_boundary_condition(grid_2d, 3)\n        \n        # Test clearing BCs\n        clear_boundary_conditions!(grid_2d)\n        default_bcs = get_boundary_conditions(grid_2d)\n        @test all(bc -> isa(bc, NeumannBC) && bc.value ≈ 0.0, default_bcs)\n    end\n    \n    @testset \"BC Validation\" begin\n        grid_1d = UnitGrid(32, (0.0, 1.0))\n        \n        # Test valid BCs\n        @test validate_bc(PeriodicBC(), grid_1d, 1) == true\n        @test validate_bc(DirichletBC(1.0), grid_1d, 1) == true\n        @test validate_bc(NeumannBC(0.5), grid_1d, 1) == true\n        @test validate_bc(RobinBC(1.0, 0.1), grid_1d, 1) == true\n        \n        # Test axis bounds checking\n        @test_throws BoundsError validate_bc(DirichletBC(1.0), grid_1d, 2)\n        @test_throws BoundsError validate_bc(DirichletBC(1.0), grid_1d, 0)\n        \n        # Test MixedBC validation\n        mixed_bc = MixedBC(DirichletBC(1.0), NeumannBC(0.0))\n        @test validate_bc(mixed_bc, grid_1d, 1) == true\n        \n        # Test TimeDependentBC validation\n        td_bc = TimeDependentBC(t -> DirichletBC(t))\n        @test validate_bc(td_bc, grid_1d, 1) == true\n        \n        # Test invalid time-dependent BC\n        invalid_td_bc = TimeDependentBC(t -> \"invalid\")\n        @test_throws ArgumentError validate_bc(invalid_td_bc, grid_1d, 1)\n    end\n    \n    # =============================================================================\n    # BOUNDARY CONDITION APPLICATION TESTS  \n    # =============================================================================\n    \n    @testset \"BC Application - Dirichlet\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        field_data = collect(1:10) * 1.0  # [1.0, 2.0, ..., 10.0]\n        \n        # Apply Dirichlet BC with value 5.0\n        bcs = [DirichletBC(5.0)]\n        set_boundary_conditions!(grid_1d, bcs)\n        \n        apply_boundary_conditions!(field_data, grid_1d)\n        \n        # Both boundaries should be set to 5.0\n        @test field_data[1] ≈ 5.0\n        @test field_data[end] ≈ 5.0\n    end\n    \n    @testset \"BC Application - Neumann\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        field_data = collect(1:10) * 1.0\n        original_interior = copy(field_data[2:end-1])\n        \n        # Apply zero Neumann BC\n        bcs = [NeumannBC(0.0)]\n        set_boundary_conditions!(grid_1d, bcs)\n        \n        apply_boundary_conditions!(field_data, grid_1d)\n        \n        # Interior points should be unchanged\n        @test field_data[2:end-1] ≈ original_interior\n        \n        # Boundary values should be adjusted for zero derivative\n        h = spacing(grid_1d)[1]\n        @test abs(field_data[1] - field_data[2]) < h  # Approximately zero derivative\n    end\n    \n    @testset \"BC Application - Robin\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        field_data = ones(10) * 2.0  # Constant field\n        \n        # Apply Robin BC: 1*f + 0.1*(df/dn) = 3.0\n        bcs = [RobinBC(1.0, 0.1, 3.0)]\n        set_boundary_conditions!(grid_1d, bcs)\n        \n        apply_boundary_conditions!(field_data, grid_1d)\n        \n        # For constant field with f=2, df/dn=0, so 1*2 + 0.1*0 = 2 ≠ 3\n        # Robin BC should adjust boundary values\n        @test field_data[1] != 2.0  # Should be modified\n        @test field_data[end] != 2.0  # Should be modified\n    end\n    \n    @testset \"BC Application - Mixed\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        field_data = collect(1:10) * 1.0\n        \n        # Different BCs on each side\n        mixed_bc = MixedBC(DirichletBC(0.0), NeumannBC(1.0))\n        bcs = [mixed_bc]\n        set_boundary_conditions!(grid_1d, bcs)\n        \n        apply_boundary_conditions!(field_data, grid_1d)\n        \n        # Lower boundary should be Dirichlet\n        @test field_data[1] ≈ 0.0\n        \n        # Upper boundary should be Neumann - hard to test without more setup\n        @test field_data[end] != 10.0  # Should be modified\n    end\n    \n    @testset \"BC Application - Time Dependent\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        field_data = ones(10)\n        \n        # Time-dependent Dirichlet BC\n        td_bc = TimeDependentBC(t -> DirichletBC(sin(t)))\n        bcs = [td_bc]\n        set_boundary_conditions!(grid_1d, bcs)\n        \n        # Apply at t=π/2\n        time = π/2\n        apply_boundary_conditions!(field_data, grid_1d; time=time)\n        \n        # Boundaries should be sin(π/2) = 1.0\n        @test field_data[1] ≈ 1.0\n        @test field_data[end] ≈ 1.0\n        \n        # Check that internal BC was updated\n        @test td_bc.current_bc[].value ≈ sin(time)\n    end\n    \n    # =============================================================================\n    # MULTI-DIMENSIONAL TESTS\n    # =============================================================================\n    \n    @testset \"2D BC Application\" begin\n        grid_2d = CartesianGrid([8, 6], [(0.0, 1.0), (0.0, 1.0)])\n        field_data = rand(8, 6)\n        original_data = copy(field_data)\n        \n        # Mixed BCs: periodic in x, Dirichlet in y\n        bcs = [PeriodicBC(), DirichletBC(2.0)]\n        set_boundary_conditions!(grid_2d, bcs)\n        \n        apply_boundary_conditions!(field_data, grid_2d)\n        \n        # X-direction (periodic): opposite boundaries should be equal\n        @test field_data[1, :] ≈ field_data[end, :]\n        \n        # Y-direction (Dirichlet): boundaries should be 2.0\n        @test all(field_data[:, 1] .≈ 2.0)\n        @test all(field_data[:, end] .≈ 2.0)\n    end\n    \n    # =============================================================================\n    # GHOST CELL TESTS\n    # =============================================================================\n    \n    @testset \"Ghost Cell Management\" begin\n        grid_1d = UnitGrid(5, (0.0, 1.0))\n        field_data = [1.0, 2.0, 3.0, 4.0, 5.0]\n        \n        # Test Dirichlet ghost values\n        dirichlet_bc = DirichletBC(10.0)\n        ghost_lower = get_ghost_values(dirichlet_bc, field_data, grid_1d, 1, lower_side)\n        @test all(ghost_lower .≈ 10.0)\n        \n        # Test Neumann ghost values\n        neumann_bc = NeumannBC(1.0)\n        h = spacing(grid_1d)[1]\n        ghost_neumann = get_ghost_values(neumann_bc, field_data, grid_1d, 1, lower_side)\n        expected_ghost = field_data[2] - 2 * h * 1.0  # Extrapolation for derivative=1.0\n        @test ghost_neumann[1] ≈ expected_ghost\n    end\n    \n    # =============================================================================\n    # UTILITY FUNCTION TESTS\n    # =============================================================================\n    \n    @testset \"Utility Functions\" begin\n        # Test BC type checking\n        @test is_homogeneous(DirichletBC(0.0)) == true\n        @test is_homogeneous(DirichletBC(1.0)) == false\n        @test is_homogeneous(NeumannBC(0.0)) == true\n        @test is_homogeneous(NeumannBC(1.0)) == false\n        @test is_homogeneous(RobinBC(1.0, 0.1, 0.0)) == true\n        @test is_homogeneous(RobinBC(1.0, 0.1, 1.0)) == false\n        \n        # Test BC string conversion\n        @test bc_string(PeriodicBC()) == \"periodic\"\n        @test bc_string(DirichletBC(1.5)) == \"dirichlet:1.5\"\n        @test bc_string(NeumannBC(-0.5)) == \"neumann:-0.5\"\n        @test bc_string(RobinBC(1.0, 0.1, 2.0)) == \"robin:1.0,0.1,2.0\"\n        \n        mixed_bc = MixedBC(DirichletBC(1.0), NeumannBC(0.0))\n        @test contains(bc_string(mixed_bc), \"mixed\")\n        \n        # Test grid utility functions\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        \n        # Default should be natural BC\n        @test has_natural_bc(grid_1d) == true\n        \n        # Set periodic BC\n        set_boundary_conditions!(grid_1d, [PeriodicBC()])\n        @test is_periodic(grid_1d, 1) == true\n        @test has_natural_bc(grid_1d) == false\n        \n        # Set non-zero Neumann\n        set_boundary_conditions!(grid_1d, [NeumannBC(1.0)])\n        @test has_natural_bc(grid_1d) == false\n    end\n    \n    # =============================================================================\n    # FUNCTION BC TESTS\n    # =============================================================================\n    \n    @testset \"Function-Based BCs\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        \n        # Test function evaluation\n        func_bc = DirichletBC((coord, t) -> coord + t)\n        value = evaluate_bc_value(func_bc.value, grid_1d, 1, 2.0)\n        # Should evaluate at boundary center (0.5) with time 2.0\n        @test value ≈ 2.5\n        \n        # Test time-only function\n        time_func_bc = DirichletBC(t -> 2*t)\n        value_time = evaluate_bc_value(time_func_bc.value, grid_1d, 1, 3.0)\n        @test value_time ≈ 6.0\n    end\n    \n    # =============================================================================\n    # INTEGRATION WITH DIFFERENTIAL OPERATORS\n    # =============================================================================\n    \n    @testset \"Operator Modification\" begin\n        grid_1d = UnitGrid(5, (0.0, 1.0))\n        \n        # Create identity matrix as test operator\n        op_matrix = Matrix{Float64}(I, 5, 5)\n        original_op = copy(op_matrix)\n        \n        # Test Dirichlet BC modification\n        modify_operator_for_bc!(op_matrix, DirichletBC(1.0), grid_1d, 1)\n        \n        # First and last rows should be identity (for boundary enforcement)\n        @test op_matrix[1, 1] ≈ 1.0\n        @test all(op_matrix[1, 2:end] .≈ 0.0)\n        @test op_matrix[5, 5] ≈ 1.0\n        @test all(op_matrix[5, 1:4] .≈ 0.0)\n        \n        # Test periodic BC modification\n        op_matrix2 = copy(original_op)\n        modify_operator_for_bc!(op_matrix2, PeriodicBC(), grid_1d, 1)\n        \n        # Should connect boundaries (simplified test)\n        @test op_matrix2 != original_op  # Should be modified\n    end\n    \n    # =============================================================================\n    # ERROR HANDLING AND EDGE CASES\n    # =============================================================================\n    \n    @testset \"Error Handling\" begin\n        grid_1d = UnitGrid(10, (0.0, 1.0))\n        \n        # Test invalid dimensions\n        @test_throws ArgumentError set_boundary_conditions!(grid_1d, [])\n        @test_throws ArgumentError set_boundary_conditions!(grid_1d, [DirichletBC(1.0), NeumannBC(0.0)])\n        \n        # Test bounds errors\n        @test_throws BoundsError get_boundary_condition(grid_1d, 0)\n        @test_throws BoundsError get_boundary_condition(grid_1d, 2)\n        \n        # Test small grids\n        small_grid = UnitGrid(2, (0.0, 1.0))\n        small_data = [1.0, 2.0]\n        \n        set_boundary_conditions!(small_grid, [DirichletBC(0.0)])\n        apply_boundary_conditions!(small_data, small_grid)  # Should not crash\n        \n        @test small_data[1] ≈ 0.0\n        @test small_data[2] ≈ 0.0\n    end\n    \n    # =============================================================================\n    # PERFORMANCE TESTS\n    # =============================================================================\n    \n    @testset \"Performance\" begin\n        # Test with larger grids to ensure reasonable performance\n        large_grid_1d = UnitGrid(1000, (0.0, 1.0))\n        large_data = rand(1000)\n        \n        bcs = [DirichletBC(1.0)]\n        set_boundary_conditions!(large_grid_1d, bcs)\n        \n        # This should complete quickly\n        @time apply_boundary_conditions!(large_data, large_grid_1d)\n        \n        @test large_data[1] ≈ 1.0\n        @test large_data[end] ≈ 1.0\n        \n        # Test 2D performance  \n        large_grid_2d = CartesianGrid([100, 100], [(0.0, 1.0), (0.0, 1.0)])\n        large_data_2d = rand(100, 100)\n        \n        bcs_2d = [PeriodicBC(), NeumannBC(0.0)]\n        set_boundary_conditions!(large_grid_2d, bcs_2d)\n        \n        @time apply_boundary_conditions!(large_data_2d, large_grid_2d)\n        \n        # Verify periodic BC\n        @test large_data_2d[1, :] ≈ large_data_2d[end, :]\n    end\n    \nend  # End main testset\n\n# =============================================================================\n# INTEGRATION TESTS WITH PDE SOLVING\n# =============================================================================\n\n@testset \"BC Integration with PDEs\" begin\n    \n    @testset \"Diffusion with Dirichlet BC\" begin\n        # Test that diffusion equation respects Dirichlet BCs\n        grid = UnitGrid(32, (0.0, 1.0))\n        \n        # Set up diffusion with fixed boundaries\n        bcs = [DirichletBC(0.0)]  # Zero at boundaries\n        set_boundary_conditions!(grid, bcs)\n        \n        # Initial condition: Gaussian in center\n        field_data = [exp(-((i-16)/5)^2) for i in 1:32]\n        \n        # Apply BCs\n        apply_boundary_conditions!(field_data, grid)\n        \n        # Boundaries should be zero\n        @test field_data[1] ≈ 0.0\n        @test field_data[end] ≈ 0.0\n        \n        # Interior should be unchanged initially  \n        @test field_data[16] > 0.5  # Peak should remain\n    end\n    \n    @testset \"Steady State with Neumann BC\" begin\n        # Test steady state solution with Neumann BCs\n        grid = UnitGrid(11, (0.0, 1.0))\n        \n        # Zero flux boundary conditions\n        bcs = [NeumannBC(0.0)]\n        set_boundary_conditions!(grid, bcs)\n        \n        # Constant field should satisfy zero flux BC\n        field_data = ones(11) * 5.0\n        original_field = copy(field_data)\n        \n        apply_boundary_conditions!(field_data, grid)\n        \n        # With zero Neumann BC, constant field should be minimally changed\n        # (some small modifications at boundaries for ghost cells)\n        @test norm(field_data - original_field) < 1.0  # Small change\n    end\n    \nend"